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Abstract 

A new lattice based scheme for numerical relativity will be presented. The scheme uses the same data as 
would be used in the Regge calculus (eg. a set of leg lengths on a simplicial lattice) but it differs signihcantly 
in the way that the held equations are computed. In the new method the standard Einstein held equations 
are applied directly to the lattice. This is done by using locally dehned Riemann normal coordinates to 
interpolate a smooth metric over local groups of cells of the lattice. Results for the time symmetric initial 
data for the Schwarzschild spacetime will be presented. It will be shown that the scheme yields second order 
accurate estimates (in the lattice spacing) for the metric and the curvature. It will also be shown that the 
Bianchi identities play an essential role in the construction of the Schwarzschild initial data. 


1. Introduction 

One of the most common techniques employed in numerical general relativity is that of hnite 
differences. In this approach the coordinate components of the metric are sampled at a series 
of well dehned points. The connection and curvatures are then calculated from well known 
hnite diherence formulae. Such formula can be easily derived in a number ways, and in 
particular, by way of a quadratic interpolation of the sampled values. Though this method 
has produced excellent results it is reasonable to explore other alternatives. 

The alternative to be developed in this paper is to use a lattice of geodesic segments to 
provide the samples of the spacetime metric. The curvature will, as with hnite diherences, be 
obtained from a quadratic interpolation of the sampled metric. This procedure is, however, 
somewhat more involved than that for hnite diherences. 

Lattice theories of general relativity such as the Regge calculus [1] are not particularly 
new. What makes our approach new is the way in which the discrete helds equations are 
constructed from the lattice variables such as the leg lengths. In the Regge calculus this is 
done by way of an action integral whereas as in our approach the held equations of general 
relativity are applied directly to the lattice. We will defer further comparison of the Regge 
calculus with our approach until after the main results of this paper have been presented. 

For any explicitly given spacetime it is straightforward, though perhaps computationally 
tedious, to construct a lattice of geodesic segments (eg. by sub-dividing the spacetime into 


a set of simplices). This lattice will inherit information about the topology and geometry 
of the parent spacetime. If the lattice is chosen correctly then it is possible to record all 
of the topological information of the parent spacetime in the lattice. The same is not 
true for the metric since the lattice contains only a hnite number of geodesic segments. 
Thus some information will be lost in going from the parent spacetime to the lattice. The 
interesting question is how much information about the parent spacetime can be recovered 
from the lattice? In particular, can the lattice be used to estimate the curvature tensor 
of the spacetime, and if so, with what accuracy? One can expect that the answer to this 
question will entail, just as in the hnite difference case, a locally quadratic interpolation of 
the sampled data. 

For the moment suppose that such an interpolation procedure exists (while putting aside 
the issue of accuracy). Then given a set of leg lengths one can compute an associated 
set of curvatures. Einstein’s equations can then be evaluated and the leg lengths adjusted 
accordingly. Clearly this amounts to solving Einstein’s equations for the parent spacetime 
discretised on a lattice. This is the basis for using lattices - they provide an alternative to 
hnite diherences as a tool for numerical relativity. 

The procedure for estimating the curvatures from the leg lengths will be based on Riemann 
normal coordinates [2-5] (later referred to as RNC). A simple algebraic dehnition of these 
coordinates is that, for a given point O, they are a local set of coordinates, such that 
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at O. 

An equivalent and perhaps more instructive geometric dehnition goes as follows. Consider 
the point O and a small neighbourhood of O. If P is any other point near to O then there 
exists a unique geodesic joining O to P. Let be the components of the unit vector to 
this geodesic at O and let s be the geodesic length from O to P. Then the Riemann normal 
coordinates of P are dehned to be = sa^. This construction fails whenever the geodesic 
joining O to P is not unique (ie. when geodesics cross). Fortunately the neighbourhood of 
O can always be chosen to be small enough so that this problem does not arise. Incidently, 
this displays the local nature of Riemann normal coordinates. They cannot be used to cover 
the whole manifold. Instead a collection of distinct Riemann normal coordinates must be 
constructed. 

From this simple dehnition a number of very useful theorems can be proved (see the Ap¬ 
pendix). In particular, for the parent spacetime, the metric is of the form 

+ 0 ( 6 ^ (IT) 
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where gfn, and R^aup are the components of the metric and Riemann tensors at the origin 
O of the RNC and and e is a typical length scale of the neighbourhood in which the RNC is 
dehned. This equation is nothing more than a Taylor series expansion of the metric around 
the origin O. Higher order approximations can be generated by continuing the series. The 
question of convergence is not very important in our context because we can always choose 
the domain of the RNC sufficiently small so that the truncation error in the above expansion 
is negligible. Note that the above expansion can be viewed as a flat space part (the constant 
term) plus a small perturbation (the curvature term). 

From this metric one can derive formulae for such things as the length of a geodesic segment 
and the angle subtended at the vertex of a geodesic triangle. The length of a geodesic 
segment with vertices i and j is given by 

L% = - ^R^aup xfx^xfx^ + C>(e®) (1.2) 

while the angle, 6k, subtended at the vertex k of the triangle with vertices i,j and k can be 
computed from 


2LikLjk cos 9k = Llk + L]k - - ^R^,aup + C>(e®) (1.3) 

These formulae can be easily applied if we are given the various g^^, R^uap and x6 of the 
parent spacetime. However, we wish to pose the inverse problem, given just the leg lengths, 
angles and the combinatoric information of a lattice, compute (or estimate) the curvature 
of the parent spacetime. This can be done by treating the above as a coupled system of 
non-linear equations. There is one such system of equations for each vertex. 

We will dehne a smooth lattice to be a lattice for which we have solved the above equations, 
(1.2) and (1.3), which in turn will be referred to as the smooth lattice equations. (Though a 
better choice might be Riemann lattice.) 

There are a number of attractive features in the smooth lattice approach that deserve spe¬ 
cial mention. First, we have made explicit use of the equivalence principle by setting the 
connection to zero at the origin of each RNC. This not only simplihes many calculations but 
it also allows us to interpret each RNC as a locally freely falling frame. This central use 
of the equivalence principle must surely be to our advantage. The second point is that the 
smooth lattice equations are very easy to formulate and solve. In contrast, the equations 
in the Regge calculus are very tedious with numerous inverse trigonometric and hyperbolic 
functions and are far from simple to setup and evaluate on a computer. Solving the Regge 
equations on a computer could easily be a very time consuming task (much more so than 
either for hnite differences or for our method). The third point is that as our method has a 
solid theoretical basis we are able to use all of the usual tools of differential geometry and 
analysis to investigate issues such as convergence and discretisation errors. Such is not the 
case with the Regge calculus. 
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Despite these attractive features one cannot apply the method naively. In any practical 

application one may need to address the following issues (amongst others). 

o Over which region should each RNC be defined? To each vertex there should be assigned a 
set of simplices from which the information is to be extracted. The size of that region (ie. 
the number of simplices) must reflect the amount of information required. It seems that 
the smallest practical region would be that composed of the ball of simplices attached to 
the nominated vertex. 

o Are there any gauge freedoms in choosing the RNC? Yes, for example the standard rota¬ 
tions, translations etc. This gange freedom can be used to force g^i, and or various xh 
to take on specihc valnes. The remaining quantities must be computed together with the 
curvatures from the smooth lattice equations. 

o Are there as many equations as unknowns? Yes, but only in two dimensions and only 
when the surface consists of triangles. In higher dimensions the set of equations is usnally 
overdetermined. This problem can be overcome by either resorting to a least sqnares so- 
Intion or by simply excluding some of the equations. Another option would be to increase 
the number of parameters, such as derivatives of the curvature, in the interpolation of the 
metric so as to prodnce a properly determined set of eqnations. 

o Is there a unique solution to the smooth lattice equations? In general, no - the smooth lat¬ 
tice equations are non-linear and so we can expect more than one solution. The members 
of a discrete family of solutions will most probably be related by discrete transformations 
such as folding one simplex over another. We would hope that there would be one excep¬ 
tional case where all the simplices do not overlap and this would be the solution we would 
choose. There is also the possibility of having a continuous family of solutions. This 
would suggest that the lattice is improperly dehned in that it lacks snfficient strncture 
(eg. too few leg lengths). This later case should not be allowed to occnr in practice. 

o Can all of the leg lengths and angles be freely chosen? If the curvatnres are specihed 
locally, as they are in our system of RNC’s, then there is qnestion as to whether or 
not there exists a single metric which has a curvature matching the specihed local values. 
This is nothing more than a question of integrability - given a curvatnre tensor does there 
exist a corresponding metric? The answer is yes provided that that curvatnre satishes the 
Bianchi identities. For our set of RNC’s this imposes (surprisingly) a set of constraints 
on the choice of the leg lengths and angles. 

Some of these issues will be explored in greater detail in later sections. 
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2. Schwarzschild initial data 


As a test of the smooth lattice method we should be able to successfully recover the time 
symmetric 3-geometry for the Schwarzschild spacetime. This metric can be written in the 
form 

where / is the proper distance measured from the throat, ri{l) is a smooth function of I and 
ig standard metric of a unit 2-sphere. 

We will present two examples of a smooth lattice discretisations of this metric. In the first 
example we will discretise the metric of the unit 2-sphere while retaining a continuous radial 
coordinate. We will first use the smooth lattice method to estimate the curvature of a unit 
2-sphere (for which we all know the exact value, 2). This estimate will then be used in a 
radial integration of the Hamiltonian constraint to complete the construction of the 3-metric. 

This is a very simple test of the method. A much stronger test will be presented in our 
second example where the full 3-metric will be discretised (though we will use the spherical 
symmetry to simplify the calculations). 

2.1. Smooth lattice 2-sphere 

The smooth lattice equations will be used in this example solely as a means to estimate the 
Riemann curvature scalar of a unit 2-sphere. This is known to have the value 2 so this may 
appear to be a somewhat un-necessary application of the method. However, if the method 
fails in this simplest of all cases then it surely will be of little use in any other situation. Once 
the curvature has been estimated we will employ it in solving the standard time symmetric 
constraint equations for the three metric. 

We will estimate R from the known 2-metric of the unit 2-sphere. To do so it will be 
necessary to choose a lattice on the 2-sphere, compute the geodesic leg lengths and hnally 
solve a set of equations (given below) for R. 

The simplest approximation to a 2-sphere is a regular tetrahedron (see Figure (1)). This can¬ 
not be expected to provide very accurate estimates for the curvature. Thus some scheme for 
rehning the lattice will be required. The scheme used in this example will be to successively 
sub-divide each triangle according to the pattern in Figure (2). 

The metric of the 2-sphere can be written as 

= d6‘^ + sin^ 6 

or as the induced metric on the surface 1 = R y‘^ + in Euclidean 3-space with the 
usual (x, y, z) coordinates. The (d, 0) and (x, y, z) coordinates are related by the usual polar 
coordinate transformations. The (x, y, z) coordinates of the four vertices of the original 
tetrahedron are easily calculated by appealing to the symmetry of the (regular) tetrahedron. 
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The coordinates of the vertices of the successive lattices are calculated in a two step process. 
First each new vertex is introduced to the centre of each old leg. This vertex is then displaced 
along the radial direction out to the unit sphere. In this way the {6, 0) coordinates of each 
vertex can be calculated. 

The leg lengths for each leg are calculated by solving the geodesic equations on the 2-sphere 
as a two point boundary value problem. At the same time we compute J ds along this 
geodesic path. This gives us the leg lengths for each geodesic segment. We then discard the 
continuum metric and turn to the smooth lattice equations to estimate R. 


2.2. Smooth lattice equations 

Consider the Riemann normal coordinate frame centred on vertex O. Suppose there are n 
triangles attached to this vertex and that the vertices, starting with O, are labelled, 0 to n. 

We are free to choose our Riemann normal coordinates such that g^v{.Xo) = diag(l, 1), {xo) = 
(0,0) and {xi) = (*,0) where -k denotes a number to be computed from the smooth lattice 
equations. This exhausts all coordinate freedoms, and thus all of the remaining x^ and 
curvature components must be computed from the given leg lengths and the smooth lattice 
equations (1.2). These equations were applied only to the legs of the triangles attached to 
O. Note that in 2-dimensions, there is only one independent curvature component, which 
we can take to be Ri 2 i 2 - 

Since there are n triangles attached to O, there will be 2n leg lengths There are also 
n + 1 vertices for which there are 2(n -|- 1) coordinates x^ to compute. However, we have 
already chosen 3 of the 2{n + 1) coordinates. Thus we have to compute 2n — 1 coordinates 
and one curvature component from the 2n leg lengths Ly. Fortunately, we have as many 
equations as unknowns. (In fact it is easy to see that this will always be true in 2-dimensions, 
provided the surface is fully triangulated.) 

The 2n equations (1.2) were solved for the x^ and R 1212 via a Newton-Raphson method. 
Starting from flat space, the iterations converged in about 3-4 iterations (though more iter¬ 
ations were required for the very coarse approximations of the original tetrahedron). 

The estimates so obtained are listed in Table 1. Since not every vertex in each approximation 
is equivalent to every other vertex (they have differing local triangulations) the method 
returns different estimates for R for each vertex. Hence in the table we have listed the best 
and worst estimates for R. One can observe that the method converges by a factor of four 
with each successive sub-division. As the leg lengths are halved with each sub-division this 
implies the error in R varies as 0{L‘^) where L is a typical length scale for leg lengths. That 
is, the smooth lattice yields 2nd-order accurate estimates for the curvature. 
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Table 1. Estimates of |i? — 2| for a unit 2-sphere 


Sub-division 

Worst estimate 

Best estimate 

Average estimate 

1 

1.19 

1.19 

1.19 

2 

2.99e-l 

4.91e-2 

1.71e-l 

3 

6.40e-2 

1.91e-3 

2.55e-2 

4 

1.54e-2 

5.67e-5 

2.97e-3 

5 

3.81e-3 

5.18e-6 

2.47e-4 


2.3. The Hamiltonian constraint 

Using the above values for R we can now proceed to the construction of the Schwarzschild 
3-nietric. Our starting point is to propose a 3-nietric in the form 

ds^ = + r]{l)^ [d6‘^ + sin^ 6 d(jP‘') 

Where / is the radial proper distance measured from the throat. There is only one non-trivial 
constraint equation, the Hamiltonian constraint, 


which must be solved for ?](/). To begin, write the metric in the 2+1 form 

where is the metric of the unit 2-sphere. Then it is a straightforward calculation to show 
that the Hamiltonian constraint is 

^^ _ 3 

dl \r] dl ) 47]“^ 2 \r] dl J 

where R is the scalar curvature of the unit 2-sphere. 

The boundary conditions at the throat, where / = 0, are chosen to be 

t] = 2 and 

The hrst condition is equivalent to setting the ADM mass to m = 1 (ie. rj = 2m). The 
second condition is required for / = 0 to be a minimal surface. 
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Note that in the above equation R can take on any value not just the discrete values found in 
the previous section. As we have already seen that the smooth lattice method gives accurate 
and convergent estimates for R, we will in the following discussions allow R to take on any 
value in the range 2 < R < 3. 

For each choice of R the initial value problem (2.3.1) was solved using a 4-th order Runge 
Kutta method starting from the throat and integrating outwards. A step length of dl = 0.2 
was used in each of the following calculations. 

The result of the integration is that we have the 3-metric in the form 

ds^ = dl^ + if {1) dVL^ (2.3.3) 

where dVt^ is the metric of the unit 2-sphere. We would like to compare this with the exact 
metric corresponding to i? = 2. One easy way to do this is to integrate (2.3.1) again but 
this time with R = 2. This will yield a metric of the form 

dS‘^ = dl^+ f{l)dQ^ (2.3.4) 

The error can then be easily computed as 

e(/,R) = -l + ?^ (2.3.5) 

r/(/) 

The results are shown in Figures (3) and (4). The hrst graph shows that for a hxed value of 
R the error rises steeply from the throat and then settles to a constant value independent of 
1. This is easy to understand, it shows that in the distant almost flat regions of the metric 
the error in approximating a smooth 2-sphere with a smooth lattice does not depend on the 
size of the 2-sphere. This graph also shows that for any hxed I the error vanishes as i? —2 
(as it must). The second plot, Figure (4), is just a series of cross-sections of the hrst plot. 
Figure (3), at specihc values of /. This graph clearly shows that the error in rj, at a hxed /, 
appears to vary linearly with the error in R. Since we have previously established that the 
error in R varies as 0{L?) we can infer that the global discretization error in rj appears to 
be C>(L2). 

It should be noted that the metrics generated for various values of R are not isometric to each 
other. To see this note that (2.3.1) is symmetric under the transformation R k‘^R, rj ^ krf 
for any k. Thus we can always scale R to R = 2. However the resulting re-scaling of rj is, 
through the boundary condition (2.3.2), equivalent to changing the ADM mass. This proves 
the assertion. The upshot of this that there is no single method to compare two metrics 
with differing values of R. The method described above is just one of many possible ways to 
compare our approximate metric against the exact metric. If one is not careful it is possible 
to generate an apparently acceptable dehnition of the error which displays, for hxed /?, an 



error that diverges as / —cx). 
form 


As an example, suppose we wrote our exact metric in the 
ds‘^ = p^(r) dO^) 


where p{r) = 1 + (m/2r). To compare this metric with our approximate metric (2.3.3) we 
could generate a transformation between the r and / coordinates by integrating, in parallel 
with the main equation (2.3.1), 


dr r 
dl 7] 


(2.3.6) 


starting from r = m/2at/ = 0. The error could then be dehned as 


The results are displayed in Figure (5) and clearly display the stated divergence. However, 
at a hxed value of I we hnd, see Figure ( 6 ), that the error vanishes as i? —> 2, ie. at each 
physical point the approximate metric converges to the exact metric. 

Note that in this example all that we have changed is the way in which the two metrics are 
compared. Why did the hrst method work so much better than the second? The answer 
lies in the choice of the transformation (2.3.6). This contains the function r] which is only 
an approximation to the exact fj. The result is that r and I are not properly aligned. What 
this means is the following. To each metric we can compute the proper distance from the 
throat to any point in question. For ds the distance to the point with coordinate I is just 1. 
For the exact metric ds the distance is /(r) = point with coordinate r. 

The above transformation (2.3.6) does not produce I = l{r). In fact, when both / and I are 
large. 



which when combined with the above (2.3.6) and dr/dl = r/fj leads to 




Thus I — I diverges (for i? 7 ^ 2) as / — 00 . This is the source of the apparent divergence 
in the metrics (for hxed R). In contrast, the original dehnition of the error (2.3.5) does not 
suffer from this problem because it compares the metric coefficients rj and f) at the same 
proper distance from their respective throats. 

An honest assessment of the above example is that even though it has given the correct 
answers it must be viewed as a very benign test of the smooth lattice method. The sole 
contribution of the smooth lattice was to aid in the computation of the scalar 2 -curvature, 
which was already known to be 2 . One could probably concoct any number of schemes which 
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spit out the magic number 2. A far better test would be to use the smooth lattice approach 
to fully discretise the 3-metric. This brings us to our second example. 


3. Fully discretised 3-metric 

Our aim in this example is to subject the smooth lattice approach to a much more stringent 
test than that used in the previous example. We will base our test on a fully discretised three 
dimensional lattice. Our plan of attack is as follows. First, we will choose the structure of 
our lattice. Coordinate and gauge conditions that reflect the desired spherical symmetries 
will then be imposed after which the smooth lattice equations will be written out in full. 
Finally we will argue that the Bianchi identities will need to be imposed so as to produce a 
set of equations from which the correct metric can be obtained. 

Consider now the construction of a lattice for a spherically symmetric space. An extreme 
example would be to choose a lattice in which the vertices have been randomly scattered 
thronghout the space. This is not only impractical but it also fails to take advantage of the 
obvious symmetries of the space. We will instead choose a lattice (see Figures (7,8)) built 
from a single tube stretching from the throat ont to the distant flat regions of the space. The 
tube has a rectangular cross-section and it is subdivided into a seqnence of cube like cells. 
Each cell stretches from one 2-sphere to the next and each snccessive pair of cells dehnes the 
region of each RNC. The fonr edges of the tnbe, not just the radial edges of the cells, will 
be required to be global radial geodesics. It is important to note that the rectangular tiles 
such as that dehned by the vertices 1,2,3 and 4, do not lie in the 2-spheres (except at the 
throat). The edges snch as (1, 2) are geodesic segments of the full 3-metric and thus cannot 
also be geodesic segments of the 2-spheres. 

Though we have only constructed one tnbe we will assnme, on the basis of spherical symme¬ 
try, that it captures all of the important geometric information of the 3-metric. Thus there 
is no need to replicate this tnbe through out the space. 

Let us now turn to the issne of coordinate and gauge conditions. This will entail arranging our 
RNC in each pair of cells, imposing restrictions on the curvature components and imposing 
the gange condition that the radial edges of the tnbe are global geodesics. 

We are free to choose an orthonormal RNC frame at the origin. Thus = diag(l, 1,1) 
at the origin. The Riemann normal coordinates for the vertices in each pair of cells were 
chosen as per Table 2. This choice can be achieved as follows. First align the ; 2 -axis with the 
radial geodesic running up the centre of the tube. The origin can then be slid up and down 
this geodesic to set the coordinates of vertices 1,2,3 and 4 to zero. Finally, by a snitable 
rotation abont the z-axis, the remaining pattern amongst the x and y coordinates can be 
achieved. 
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Table 2. Coordinates for the vertices in Figure 8 


Vertex Vertex Vertex 

1 { a , a , b ) 1 ( a, a, 0) I"*" ( ^ 

2~ {a , —a , b ) 2 {a, —a, 0) 2’’’ ( a'^, —a'*’, b^) 

3~ (—a ,—a ,b) 3 {—a,—a, 0) 3'*’ (—a'*’,—a'*’, b^) 

4~ {—a , a , b ) 4 {—a, a, 0) 4’’’ (—a'*', a'*', &'*') 


The spherical symmetry of the space must impose some restrictions on the curvatnre com¬ 
ponents in the Riemann normal coordinates. Recall that the general form of a RNC metric 

is 

~ 9jJ,v ^ -\- 0^6 ) 

The only parameters which we can play with are the R^avp- Thus to respect the required 
symmetries we must restrict these parameters accordingly. To do this we revert for the 
moment to a generic spherically symmetric metric 

ds^ = dr‘^ + rj(r)‘^ {dO^ + sin^ 9 d(j?) 


where ? 7 (r) is any smooth fnnction and r is the radial proper distance measnred from the 
throat (this was written as / in the previous sections). In these coordinates there are only 
two non-trivial frame components, 


R - . = R . - = - 1 ^ 
^fOfO ^f4>r(j) ^ ^^2 



(3.1) 


where, for example, R§^§^ = Thus we expect only two non-trivial curvature 

components in the RNC frame. If we align the r axis with the axis then it is easy to see 
that the corresponding non-trivial RNC components must be Rxyxy and Rxzxz = Ryzyz- 
Thus, the RNC metric can be reduced to just 


ds^ = dx^ -|- dy"^ -|- dz^ - Rx{xdy — ydx)"^ - Rz{xdz — zdx)"^ - Rz{ydz — zdy)" 


where Rx = Rxyxy and Rz = Rxzxz = Ryzyz- Incidently, for any pair of vectors and 
we have 


R^^aup = Rxiv^uV - + Rz{v^u^ - vVf + Rz(vV - v^uyf 
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Using this metric and the above choice of coordinates we obtain, from the smooth lattice 
equations, 


0 = 3L2 - I2a^ + (3.2) 

0 = 3(L“)2 - 12(a“)2 + 4i?^(a“)^ + AR,{a~b~f (3.3) 

0 = - I2{a^f + 4i?^(a^)^ + AR^{a^h'^f (3.4) 

0 = 3(d“)2 - 6(a - oTf - 3{b~f + 2R^{ab~f (3.5) 

0 = 3{d^f - 6(a - a^f - 3(6^)^ + 2R^{ab'^f (3.6) 

0 = cos/5 + 4i?2,a^ (3.7) 

3Ld cos a = 6a{a — a )—ARxa^{a — a ) — 2Rz{ab )^ 

3Ld^ cos a'^ = 6a(a — a’*’) — AR^a^la — a’*’) — 2Rz{ab'^)‘^ 


We have previously stated that the radial edges of our tube must be global geodesics. Thus 
there can be no kink in the incoming and outgoing edges at vertices 1,2,3 and 4. This implies 
that 

CX CX — TT 

which using the above leads to 

0 = d^ (6a{a — a ) — AR^a^ia — a ) — 2Rg{ab )^') 

_ / X (3.8) 

+ d f6a(a — a'*’) — AR^a^^a — a’*’) — 2Rg{ab^)‘^] 

There are seven quantities that we must compute a, a’*’, a ,6'*’, 6 ,Rx and Rz- We thus 
need seven equations which we take to be equations (3.2-3.7) and (3.8). Though this may 
seem like a well dehned system with seven equations for seven unknowns there is a serious 
problem. Suppose we were to freely choose all of the d’’’, d~, LX, L, L~ and cos/5. We could 
then solve the lattice equations and thus obtain the curvatures in each RNC frame along the 
tube. Since this is a well dehned set of equations we see that this is equivalent to having an 
arbitrarily specihed set of curvatures along the tube. However, in the continuum metric the 
two curvatures are derived from one function and therefore they cannot be freely specihed. 
Indeed from (3.1) we see that 


0 = 



(3.9) 


where R^ = Rq^q^ and R^ = Rf.§f§- This is nothing other than the standard Bianchi identity 
for the continuum metric. We can expect that a similar but discrete version of this equation 
must also exist in our smooth lattice model. 
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3.1. Bianchi Identities 


It is well known [6] that a necessary and sufficient condition for the existence of a metric 
g^u{x) from which a given curvature tensor can be derived is just the Bianchi identities 

0 ~ Rfj,uap-,p RpLvPp\a Rpvpa\P 

In each of our RNC frames this can be written as 

0 RpuaP,p RpvPp,a ^pvpa,P 

since = 0{e^) throughout each RNC (in the conformal metric, see the Appendix). 

Unfortunately we cannot make much use of this equation as it stands. The problem is that 
the inversion of the smooth lattice equations returns the curvatures but not their derivatives. 
One might object by arguing that for the RNC metric we can evaluate the curvature and 
its derivatives to any order. This may be so but such calculations cannot be expected to be 
accurate estimates of those quantities for the smooth metric (that the lattice is interpolating). 
The same situation arises in the much simpler case of a piecewise quadratic interpolation of 
a function of one variable. The third derivatives of each interpolant is zero yet it is most 
unlikely that the function has a zero third derivative at each point. The normal practice is 
to use the second derivative of the interpolant as an estimate at just one point, usually at 
the central point. Estimates of higher derivatives can then be obtained by interpolation of 
this derived data. 

This same philosophy seems appropriate for our lattice. The estimates for the curvatures 
are valid at just one point, the origin of the RNC. Derivatives of the curvatures should be 
estimated by differentiation of a local interpolation of the point estimates of the curvatures. 
There is one slight subtlety here - the curvatures are dehned in different RNC frames. This 
will require some coordinate transformations to bring the local curvatures into one common 
frame before the interpolation and differentiation can be performed. This is actually not 
very difficult. Suppose for example that the origins of two neighbouring frames are O and 
O'. Suppose that the associated metric and curvatures are gpu,Rpvap and g^'p', R^.'pia'p' 
respectively. We seek the representation of R^'p'a'P' the RNC at O. To do this start with 
the set of coordinate basis vectors at O' and parallel transport them to O. They will be 
related to the basis vectors at O by some transformation matrix p. Then the transformed 
value of the curvature will be Rf^'p/^'p'R^ pR'^ p- For this calculation we can, to 
leading order in e, evaluate A^ p solely from the flat space parts of the metrics gfj,p and g^'p/. 

Though the procedure just described may be suitable for a generic lattice there is a simpler 
approach for our symmetric lattice. The idea is to use the equivalent integral representation 
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of the Bianchi identities. In three dimensions this happens to be 


0 


IdM 


{Rlivaf5,p Rpvpa,}^ dx dx^dx^ 

Rpuaf3 dx'^dxd 


where M is an arbitrary three dimensional region with boundary dM. To leading order in 
e the limits of integration can be set as if M was everywhere flat. This later form is much 
easy to apply since it does not require estimates of the derivatives of the curvatures. For our 
lattice we will choose M to be the pair of cells that dehne the typical RNC frame. Though 
this consists of two cells and thus ten boundary faces we can treat them, to leading order in 
e, as just one larger cell (ie. a linear cell, with six faces, that expands from L~ to L'*’ over a 
length d~ + d^). The above equation can then be written as 


0 



Rpvaf] dx dx^ I 
Rpvap dx dx^ -\- I 

Js+ 

Rpvap dx dx^ -\- I 

I 0 + 


Rpva.j3 dx dx^ 
RpvajS dx dx^ 
Rpvap dx dx^ 


where S~,S^ are the faces of the cells on the +x and —x axes respectively (with similar 
dehnitions for the remaining four faces). We will now compute the curvatures on each of 
these faces in this RNC frame. This information can be obtained by an interpolation and 
coordinate transformation of the curvatures from neighbouring cells. 

At this point, it is appropriate to represent all of the curvature tensors in bivector form. Let 
Upu, and W^p be an orthonormal set of normalised bivectors (1 = U^pU^’^, 0 = U^pV^‘' 
etc.) associated with the xy,xz and yz planes respectively. Then the curvature tensor at O 

is 

RpvaP ~ RxUp.vU + Rz^pv^afi T Rz^pv^aP 
while at we would have, in the coordinates for O^, 


RpuaP = RxKpKp + R>pXp + RtK-W^p 


However, it is clear that, to first order in e, the transformation matrix p for these two 
frames is just the identity matrix. Thus the linear interpolant for the curvature, along the 
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axis between O and is 


Rfii’apid) 


~ (d'*') ^fJ-vad 


where d is the proper distance measured along the 2 : axis from z = 0. For the face S'^ we 
have d = d'^. 

For the interpolation between O and O^, the procedure is much the same with the exception 
that the transformation matrix is now a rotation. Indeed it is easy to see that the 
bivectors in the two frames are related by 


= Uf,u cos Ap + W^i, sin Ap 
— V 

COS Ap-U^u sin Ap 

where 2sin(Ap/2) = (L'*’ — L~)/(d^ + d~). The interpolation can then be obtained by 
generalising these equations to be active transformations of the original bivectors U, V and 
W into a moving set of bivectors, namely 


cos p + sin p 

v^Ap) = 

^^v{p) = Wpi. cosp -U^usinp 


leading to 


VaMd) = RxUlAp)Kp{p) + Rzv;i,{p)v^p{p) + R,w;,{p)w 

The interpolation to the face S'^ can now be obtained by setting p = Ap/2. Note that we 
can use the approximation Ap/2 ^ sin(Ap/2) since the four radial geodesic edges of the 
tube must be close in order for the smooth lattice to be a good approximation to the smooth 
continuum metric. 

Similar interpolations can be constructed for the remaining faces. These can then be sub¬ 
stituted into the Bianchi identity with the result, assuming each is constant on each 
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face, 


0 ~ {^x 

+ R,At (+Ap/2) + H'^-,(+Ap/2)) - R^A' (\'«,(-Ap/2) + H'^ (-Ap/2)) 


where A^, and are the areas of and S'^ respectively. The areas are easily 

seen to be 



z^x — ^x 


- (l+ + L ) (d+ + d ) 


Contracting the above equation with Ufj^p and using Ap = (L’*’ 
to 


0 = {L^Rpcf - + T") 


-L )/ (d’*’ + d ) leads directly 

- L-^ Rz 


This is our discrete version of the Bianchi identity. It clearly resembles the previous Bianchi 
identity (3.9) when expressed in hnite difference form. Note that contractions with the other 
bivectors, V and W, leads only to trivial equations. 


3.2. The smooth lattice equations 

Its been a rather long road but we can now write down all of the smooth lattice equations, 
namely. 


0 =120^ - - 3L2 (3.2.1) 

0 =12(a“)2 - 4R^{a~)^ - 4R,{a~b~)‘^ - 3{L~f (3.2.2) 

0 =12(a'')2 - 4Rp,{ay - 4i?^(aV)2 - 3{L^f (3.2.3) 

0 =6(a - a~f + 3{b~f - 2R^{ab~f - 3{d~f (3.2.4) 

0 =6(a - a^)2 + 3(6'^)2 - 2R^{ab'^f - 3{d'^)^ (3.2.5) 

0 =4i?a,a^ + cos/5 (3.2.6) 

0 = d^ (6a{a — a ) — 4Rxa^{a — a ) — 2Rz{ab )^j 

+ d ^6a(a — a’*’) — 4i?a;a^(a — a”^) — 2i?2(a&’'’)^j (3.2.7) 

0 + 2Rz (3.2.8) 

0 = - {L^R^y - [l-^ + L“) - L“) Rz (3.2.9) 


The second last equation is just the Hamiltonian constraint ^^^R = 0. 
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3.3. Solution strategy 

The smooth lattice was constructed by starting from the throat and successively adding cells 
to the end of the tube. There are two parts to this scheme, first the choice of initial data at 
the throat and second, the repeated solution of the smooth lattice equations. 

Consider first the generic step in which one extra cell is added to the tube. Suppose we 
are given values for R~ and Rx (from previous calculations). Then we can solve the nine 
equations (3.2.1-3.2.9) for the nine quantities a, , a~, 6'*’, b~, Rz, R^, and cos (3. Then 
the R~ and Rx of the next RNC frame can be taken to be Rx and R^ of the current RNC 
frame (the transformation matrix between the pair of frames is just the identity matrix, to 
first order in e). This process can then be repeated as many times as might be required. 

All of the equations were solved using a Newton-Raphson method with initial guesses chosen 
to be flat space. That is, given d~, d'*’, L~ and L, all of the parameters were estimated using 
Rx = Rz = 0. This is very easy to do, leading to 


cos /5 = 0 



At the throat we need to make a few minor alterations to the above. Our boundary condition 
is that the throat must be a minimal surface. This constraint can be imposed via a standard 
reflection symmetric condition. That is, for the RNC centred on a point on the throat we 
demand that cells on either side of the throat be mirror images of each other. This leads to 

0 = a'^-a" (3.2.2') 

0 = 6+ + 6“ (3.2.4') 

and 

L+ = L~ 

= d 

This however leads immediately to a problem with the Bianchi identity at the throat in that 
it cannot be solved for a unique value for R^. There are two solutions to this problem. First, 
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one can set 


0 = R^-Rx (3.3.1) 

This follows from the observation that in the continuum limit, as d'^ ^ 0 and d~ ^ 0, the 
Bianchi identity can be reduced to 


0 = 



Rz 


which, with 0 = dL/dl at the throat, leads to 0 = dRx/dl and hence the estimate R^ = Rx- 

The second solution is to combine the above the differential equation with the constraint 
0 = Rx + 2Rz. This leads directly to 


0 ^ I iL^R.) 

and thus 

0 = {L^Rxf - {L^Rx) (3.3.2) 

Though this is valid along the whole tube we only apply it at the throat so as not to 
compromise the independence of the smooth lattice equations (we want the method to stand 
on its own and not to be supported or dragged along by results from the continuum). 

At the throat we replaced the equations (3.2.2,3.2.4) and (3.2.9) with the equations (3.2.2', 
3.2.4') and (3.3.2). The nine equations were then solved for the same nine quantities using 
as initial guesses 


cos /5 = 0 
L~ = L'^ = L 
a = a = a = -L 

—b = = d^ 

Rx = Rz = 0 

Using these initial guesses, the Newton-Raphson method for both the generic and initial 
problems converged in less then five iterations. 

Now let us return to the first point raised at the beginning of this section, namely, the choice 
of initial data L, and Rx on the throat. Collectively this corresponds to setting the ADM 
mass and the fraction of a 2-sphere covered by each tile (such as (1,2, 3,4) in Figure (8)). 
These could be chosen randomly but instead we chose to compute them from the known 
3-metric so that we could explore the convergence properties of the smooth lattice solutions. 
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Indeed the easiest way to compare the smooth lattice solutions against the exact 3-metric 
is to mimic the construction of the lattice using the exact 3-metric. Thus we construct a 
series of concentric 2-spheres and four radial geodesics. The vertices of the lattice arise as 
the intersections of the radial geodesics with the 2-spheres. All of the leg lengths, angles and 
curvatures can then be computed and compared with those from the smooth lattice. 

The exact metric in isotropic coordinates is 

= p(r)^ (dr^ -f r^(dd^ -f sin^ 6 d4>^)) (3.3.3) 

where p{r) = 1 -|- m/ (2r). For all of our calculations we chose m = 1. 

We can chose the coordinates of the vertices according to the following table. 


Table 3. Coordinates for the vertices of the exact lattice. 


Vertex 

(r, d, 0) 

Vertex 

(r, d, 0) 

Vertex 

(r, d, 0) 

1" 

1 

> 

o 

1 

(r,Ad, 0) 

1+ 

(r+,Ad, 0) 

2" 

(r ,Ad, 7r/2) 

2 

(r,Ad, 7r/2) 

2+ 

(r+,Ad, 7r/2) 

3" 

(r",Ad, tt) 

3 

(r. Ad, tt) 

3+ 

(r^Ad, tt) 

4~ 

(r~, Ad, 37r/2) 

4 

(r. Ad, 37r/2) 

4+ 

(r"^. Ad, 37r/2) 


The size of the patch on the initial 2-sphere is determined by the value of Ad. We chose 
Ad = e(7r/40) where e is a scale parameter in the range 1/256 < e < 2. The effective number 
of tiles on each 2-sphere can be estimated as N = 47r/(Ad)^ which for the chosen range, 
corresponds to 509 < N < 134 x 10® tiles. 

The values of r'*’ are obtained by solving, via a Newton-Raphson method, the equation 


+ 



= F(r^) - F(r) (3.3.4) 

where F{u) = u + m\og{2u/m) —m?/{Au). On the throat we chose r = m/2 and r~ so that 
d'^ = d~ = F{r)-F{r~). 

The exact leg lengths, which we will write as L* and L^, require a little bit of work to 
compute. They are the lengths of the geodesic segments joining the various vertices. This 
entails the solution of a two point boundary value problem 

d^x^ . X dx® dx^ 

(3.3.5) 

xt^{0) = , a;^(l) = 
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for the geodesics x^{X) for which we used a simple shooting method built on a 4-th order 
Runge-Kutta routine. The lengths were then computed by evaluating f ds along the geodesic. 
This was done by adding one extra differential equation to the Runge-Kutta routine, namely, 


dL 

dX 



dx°‘ dx^ 
dX dX 


( 1 / 2 ) 


(3.3.6) 


In solving the above two point boundary value problem we used isotropic x, y, z coordinates 
so as to avoid the coordinate singularity at 6 ^ = 0. Though we did not need to integrate 
through this singularity its presence degrades the quality of our estimates for the geodesics 
and the leg lengths. 

The frame components of the curvature are easily computed from the above metric (3.3.3) 
with the result that 


(-R®)* '^a/3d/3 

(Rz)* = Rfafa 


2m 

r^p{r)^ 

R . 

f/3r/3 r^p(r)® 


(3.3.7) 

(3.3.8) 


where p(r) = 1 -|- m/ ( 2 r). 

The above equations (3.3.4-3.3.7) were used to compute initial values for L*, L'l and {Rx)~k- 
These were then assigned to L, and Rx- The radial lengths leg lengths d^ were chosen to 
be d^ = yL where 7 is a hxed constant. This choice ensures that the cells do not become 
too fat nor too elongated as successive cells are added to the tube. This scheme provides 
for an adjustable step length, short steps near the throat where the curvatures are large but 
longer steps away from the throat where the curvatures are weak. In all of our calculations 
we set 7 = 0.1 which we arrived at by direct experimentation (this point will be discussed 
further in the following section). 


3.4. Results 

The smooth lattice equations (3.2.1-3.2.9) were repeatedly solved to generate data out to 
/ = 50. As each new cell of the lattice was constructed, the equation (3.3.4) was solved so 
that the smooth lattice quantities L, Rx and Rz could be compared with the exact values 
L*, {Rx)-k and (R^)*. 

The errors between the exact and smooth lattices depend upon the distance along the edge 
of the tube measured from the throat I = ^d^ and upon the choice of A0. The fractional 
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errors were computed as 


EL{1, e) — —1 + — 

J-i-k 

\Ez)* 

where e = ( 4 O/ 7 r)A 0 is the scale parameter. These errors are displayed in Figure (9) and 
clearly show that the errors vanish, at hxed /, as A0 —0. However, for a hxed A0, the errors 
appear to grow without bound for increasing values of 1. This is much like the behaviour 
we encountered in section (2.3) (though much less pronounced). We can infer from that 
section that the divergence in the errors in this example may also be due to an ambiguity 
in correlating the radial positions between the numerical and exact metrics. For our smooth 
lattice there are at least two ways to measure the radial proper distance. One can measure 
it along the outer edges of the tube, as we have done, or along the central geodesic (ie. along 
the axis). It is an arbitrary choice as to which of these is to be used to dehne a radial 
position in the exact metric. Both choices agree at the throat but progressively diverge as 
we move away from the throat (because the tube is always increasing in cross section). 

We have also estimated the rate of convergence by plotting the errors as function of e for 
a series of hxed values of /. These are displayed in Figure (10) and clearly show that the 
convergence of each of the quantities, L,Rx,Rz to their continuum values is quadratic in 
the lattice spacing. The wiggles in the lower left corner are probably due to round-off errors 
becoming signihcant. 

3.5. Variations 

Three variants of the above smooth lattice method were tried. First we repeated the above 
calculations but this time using equation (3.3.1) instead of (3.3.2) as our boundary condition 
on the throat. The results are plotted in Figures (11,12) and the only signihcant diherence 
is the appearance of an instability in the errors of the curvatures. The amplitude and 
wavelength of the oscillations in the curvatures were both found to decrease with decreasing 
7 . For 7 2 the errors peaked at about 40 (very large indeed!) with only about two cycles 

occurring. By trial and error we judged 7 = 0.1 to give the most pleasing results (to the 
eye). Note also that in the regions where the relative errors ER^ are large and oscillatory, 
the absolute errors Rx — (Rx)* are of order 10“^. Thus this instability is not a cause for 
great concern. Indeed it appears to have little ehect on the errors in the leg lengths (compare 
Figures 9 and 11). 

For our second variation we took advantage of the integrated form of the Bianchi identity 
(3.3.2) to replace the smooth lattice equation (3.2.9) with the equation 

k = L^Rx (3.2.9') 
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with k a constant computed on the throat from the exact metric. The nine equations 
(3.2.1 — 3.2.8) and (3.2.9') were solved for the nine quantities a, a’*’, a~, 6’’’, h~, Rz, Rx, (-h'*')^ 
and cos/5. The results for this variation were much the same as for the original method. 
However, this variation showed no signs of instability in the errors in the curvatures. 

Our third and hnal variation was designed so as to remove the need for the Bianchi identities. 
Recall that the Bianchi identity was required to constrain the two degrees of freedom in the 
curvatures Rx and Rz so that they were derivable from one metric function. We know that 
in the exact metric each tile such as (1,2, 3,4), when viewed as a tile on the 2-sphere, is a 
scaled image of the tile on the throat. It is not unreasonable to expect that when the tiles 
are built with respect to the full 3-metric a similar though approximate scaling could be used 
and that the associated errors would be negligible for small tiles. To investigate this opinion 
we altered the structure of the lattice to include a diagonal brace joining vertices 1 to 3. 
Similar braces were also included for the pairs of vertices 1~,3~ and l'*’,3'^. The lengths of 
these diagonals will be denoted by m, m~ and m’*’ respectively. The cos /3’s are discarded in 
favour of these diagonals. It is a straightforward calculation, from the main smooth lattice 
equations (1.2), to show that 

0 = Sa‘^-m^ (3.5.1) 

The condition that each tile is a scaled version of the tile on the throat can be written as 

0 = (3.5.2) 

where k' is constant along the tube. There are two apparently reasonable choices for the 
constant kh It could be calculated either from the exact 3-metric on the throat or from 
the distant regions, far from the throat, where the metric is flat. In fact both choices are 
unacceptable. The former choice, setting k' on the throat, cannot be used because it does 
not allow the lattice (with our assumed coordinates, see Table (2)) to become asymptoticly 
flat. For smaller and smaller sizes of the initial tile on the throat we can expect that the 
asymptotic curvature might also get smaller and smaller yet it seems unlikely that this could 
save the day. In the second choice, which requires {k')‘^ = 1/2, we hnd by direct inspection 
of the smooth lattice equations that 0 = Rx = Rz everywhere (which when coupled with 
the boundary condition on the throat then forces L~ = L = along the tube). Thus this 
second choice certainly cannot be used. We thus opted for the hrst choice (but with few 
expectations). 

In a typical RNC we computed, for this variation, a, a'^, a~,b'^,b~, Rz, Rx, and 

This required nine equations which we took to be the smooth lattice equations (3.2.1-3.2.9) 
with (3.2.6) (which dehnes cos/5) and (3.2.9) (the Bianchi identity) replaced by (3.5.1) and 
(3.5.2). The results are displayed in Figure (13) and are as expected - the method does not 
converge in any sense. It appears that forcing the tiles of the smooth lattice to observe the 
same symmetries as the tiles of the 2-spheres is incorrect. We conclude that the Bianchi 
identities cannot be ignored - they are essential in providing the correct constraint on the 
L’s and cos/5’s (or the L’s and m’s). 
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3.6. Are the Bianchi identities essential? 


Much has been made of the Bianchi identities and the role that they play in the smooth 
lattice equations. Yet does it not seem odd that in the smooth lattice approach we find our¬ 
selves computing one of the curvatures by integrating a differential equation yet all standard 
formulas give the curvatures as explicit functions? As another explanation as to why this is 
so consider the frame components of the curvature of the exact metric. These are 

1 cPri 
f3 dr‘^ 



{Rx)* = 





2 


If we now choose the leg lengths according to L{r) = r]{r)A6 then we get 


(R) 

L dr^ 


{Rx). = 


1 

I2 



Now suppose that all we are given are the discrete leg-lengths and that we wish to extract 
the curvatures from the discrete lattice. We can expect that in some continuum limit we 
would recover the above formulae. But where in this process would the number A6 come 
into the game? How is it extracted from the discrete lattice (prior to the continuum limit) 
and how do we know that it is a constant along the tube? The answer must lie partly in 
the constraint that the radial edges be global geodesics and partly in the use of the Bianchi 
identities. In fact it is easy to show (see section (4) below) that, in the continuum limit, 
equations (3.2.1-3.2.7), can be reduced to just 


Rz 


1 d^L 
L dr^ 


while from the remaining equations (3.2.8-3.2.9) we can obtain 


R 


X — 


1 

Z2 




where k" is a constant along the tube. That this constant equals (A9)‘^ can be seen by 
evaluating this expression for large r where the curvatures vanish. Thus using a metric such 
as (3.3.3) introduces global information that is simply not accessible from the purely local 
calculations used in one RNC of the smooth lattice. 
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This again shows the essential nature of the Bianchi identities in producing the correct 
estimates for the curvatures. 


4. The Continuum Limit 


Though the numerical results presented in the previous sections provide strong evidence 
that the smooth lattice equations are a convergent approximation to the correct differential 
equations it would be nice to establish this by direct analytical methods. 

The most direct route to the continuum is to introduce a scaling for all the leg lengths such 
as L —eL, ed^ etc., and to then express the solutions for a,d^, ■ ■ ■ as power series 

in e. The terms in the power series can be obtained by substitution into the smooth lattice 
equations. From equations (3.2.1-3.2.5) we obtain, assuming Rx and Rz are known. 


a = -^ (l + + 0(6-)) 

^ (l + + gfl. ((dp" - \(L - Ly'^ + ©(s'*)) 

V = F (^1 + - 1(L - L-)^ + oiy'j 

{by = (dy - ^{L - ly+oy) 

{b-f = {d-f-^{L-Ly+o{y 

The expansions have been carried as far as is necessary to obtain the continuum limit. These 
can be substituted into equation (3.2.7) which can then be reduced to a differential equation 
by employing the usual Taylor series expansions 


L~ 






+ ■■■ 

+ ■■■ 


where I is the proper distance measured from the throat. The result, after cancelling a 
common factor of + d~)L/2, is 

0 = + RzL^ + O(e^) + 0{e){d^ - dT) 

The hrst error term O(e^) arises from the curvature terms in the smooth lattice equations 
while the second error term 0{e){d^ — dT) arises from the Taylor series expansion in L. 
Clearly in the limit when e —0 we recover the correct differential equation. We can also see 
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that when the successive are chosen to be a smooth function of L then — d = 0{e^) 
and thus both error terms are 0{e^). 

Equations (3.2.8-3.2.9) can be treated in a similar fashion leading to 

° " Ts + (^) 

In both of the above equations the errors terms are Cl(e^) relative to the main terms. Thus 
we can infer that the relative errors EL{l,e) and ERx{l,e) should both be This is 

exactly what we observed in our numerical studies. 

We have not used equation (3.2.6) since it can be viewed as a dehnition of cos/5. 


5. The Regge calculus 

The Regge calculus [1] is another lattice theory of General relativity. Though it employs 
a lattice it differs signihcantly from our method in the way in which the held equations 
are computed from the lattice variables. In the Regge calculus the metric interior to each 
pair of cells is deemed, a-priori, to be hat. This precludes the use of standard methods 
for computing the curvatures (such as pointwise diherentiation of a metric). Instead one is 
forced to look at integral quantities, in particular the action integral 

I = [ d'^x 


The curvature for the piecewise hat metric can be interpreted in the sense of distributions 
as a series of Dirac delta functions on the 2-dimensional sub-spaces of the lattice (usually 
the triangles of a simplicial lattice). The action integral can then be evaluated as 


1 = 25 ^ OiAi 

i 


where the sum includes all of the triangles and Ai and 6i are the areas and defects on these 
triangles. The Regge held equations are then obtained by extremising the action with respect 
to the leg lengths, leading to 


0 




dLj 


There is one such equation for each leg Lj. 

Though this derivation mimics that used in Einstein’s theory one cannot state, by this 
fact alone, that the Regge calculus is a consistent discretisation of Einstein’s equations. By 
consistent we mean that as the Regge lattice is successively rehned the solutions of the Regge 
equations converge to solutions of Einstein’s equations. Proving this is extremely difficult 
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primarily because the metric is not differentiable as a point function. Another complication 
is that any metric can be approximated by many different sequences of lattices yet one would 
want to develop theorems that account for this vast array of lattices. One hrm resnlt has 
however been proven, [7], that the Regge and Hilbert actions converge when evalnated over 
a hxed region of the limit smooth space. The proof is very detailed and far from trivial. 
The state of play is (i) that any smooth metric can be approximated, to any accnracy, by a 
simplicial metric (snch metrics need not be solutions of their eqnations), (ii) the Regge and 
Hilbert actions converge when evalnated on a fixed region of the limit smooth space (again, 
for any smooth metric) bnt (iii) their remains an open qnestion as to whether or not the 
limit metric of a seqnence of solntions of Regge’s eqnations is also a solntion of Einstein’s 
eqnations. 

Of conrse the later qnestion can be tackled via direct nnmerical experiments. In all cases 
reported to date it appears that the Regge solntions do appear to converge to solntions 
of Einstein’s eqnations. However Brewin [8] has argned that this may be because the ex¬ 
periments have been condncted on highly symmetric spacetimes with the Regge lattices 
constrained by those symmetries. It could be that the symmetries are so strong as to bring 
abont a concordance between the two theories. Brewin argued that in less symmetric space- 
times the distinction between the two theories may become apparent. He went further to 
show that even for the symmetric spacetimes, snch as Schwarzschild, a lattice that was bnilt 
pnrely from 4-simphces and the known exact metric, did not appear to have the appropriate 
trnncation errors when evalnated on the Regge eqnations. This and other examples were 
used to argue that the Regge equations are fnndamentally not a consistent discretisation of 
Einstein’s eqnations. This is not a view shared by the rest of the Regge calculus community. 

In contrast there can be little argnment regarding the correct held eqnations for the smooth 
lattice approach. They are just the usual Einstein eqnations. 

Of particular note for the cnrrent discussion are the calculations of Wong [9]. He used 
the Regge calculus to solve the time symmetric initial valne eqnations for a Schwarzschild 
spacetime. He chose two methods, one in which the 2-spheres were triangnlated by icosahedra 
and a second in which he used inhnitesimal tiles to generate a rectangnlar tnbe similar to 
onr second main example. In this second method Wong proceeded to the limit of vanishingly 
small cross section in the tnbe which he called the continnnm method. However, in both 
cases Wong chose to retain a fnlly discrete radial snb-division. The results for both methods 
were very good with the typical errors in the range of 1-10% for the icosahedral method and 
0.001-1% for the continnnm method (for qnantities snch as volnmes and distances from the 
throat). 

It is tempting to compare these results with our own. In the icosahedral method there are 
20 triangles per 2-sphere. This lies between the second and third snb-divisions of onr first 
example and to a choice of e ~ 10 in onr second example. In both cases we hnd the error in 
the leg lengths to be approximately 20% while the errors in the curvatures are approximately 
50% for the first example and 100% (yikes!) for the second example. These are larger than 
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Wong’s results which is slightly surprising since htting a smooth metric to a lattice could be 
expected to be more accurate than htting a piecewise hat metric to the same lattice. 

Wong’s continuum can only be compared with our tube method. We can see from Figure 
10 that our errors are around 10“^% whereas Wong’s best results are around 10“^%. This 
comparison is a bit unfair on both methods. In all of our calculations the edges of our cells 
were scaled via the parameter e whereas Wong chose to retain a discrete radial sub-division. 
This gives an advantage to our method. However where Wong obtains an advantage over our 
method is in his extensive use of the target continuum metric. Indeed he writes the metric 
in the form 

-I- r]{r)^ {dO'^ + sin^ 6 d(jP‘') 

(with a slight change in replacing his p(r) with our ? 7 (r)) and estimates all of the leg lengths 
in the Regge lattice in terms of the unknown function ri{r) and freely chosen quantities such 
as Ad and A0. In regulating the choice of leg lengths via this metric he has avoided the 
problem of having too many degrees of freedom in the lattice. We had the same problem 
but we resolved it by way of the Bianchi identities. 

It is unlikely that such a metric ansatz will be available in less symmetric spaces forcing 
one to look at fully discrete ways in which to impose the desired symmetries. Yet an overly 
enthusiastic imposition of the symmetries on the space of possible lattices can restrict that 
space to a set of measure zero. For example, suppose we set out to use the Regge calculus 
with tubes much like that used in our second example. We could take four such tubes and 
join them together along one common radial edge. In this model we would not need the 
cos/5’s. All of the defects can then be easily calculated. In fact we hud that they are all 
non-positive. The only solution of Regge’s equations is therefore flat space. This can even 
be seen without doing any calculations - try htting the four tubes together in the distant hat 
regions where the defects should all be close to zero while dL/dr is constant but non-zero. 
Thus some of the symmetries amongst the legs must be relaxed. But this then introduces 
extra degrees of freedom. How are they to be constrained? This is where Wong’s ansatz 
saves the day. It produces sufficient variation in the leg lengths of the four tubes so as 
to obtain meaningful solutions of the Regge equations while retaining the desired spherical 
symmetries. 

It is hard to compare Wong’s continuum method and our second example on a level play¬ 
ing held as Wong’s method obtains signihcant advantage in using a method not generally 
available in a fully discrete implementation of the Regge calculus. 

The smooth lattice method could have been applied directly to Wong’s icosahedral method. 
We chose not do so because it seemed unlikely to add anything new to our investigations. 

It is important to note that the smooth lattice method uses the same lattice structures 
as used in the Regge calculus. Thus it inherits all of the attractive features of the Regge 
calculus, in particular the natural division between the topological and metric information 
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and that the metric is recorded in pure scalar form (eg. the leg lengths). But more impor¬ 
tantly, the smooth lattice method has the additional advantage that it extracts a smooth 
metric from the lattice and thus that the Einstein equations may be applied directly to the 
lattice. There is also a practical advantage of the smooth lattice method - the equations 
are much easier to construct and evaluate than those of the Regge calculus (witness some 
of the unwieldy calculations of defect angles with their attendant inverse trigonometric and 
hyperbolic functions, its a mess!). The smooth lattice equations can reasonably be expected 
to be much quicker to evaluate and to solve on a computer than the Regge equations. 


6. Discussion 

On the limited evidence provided here it is hard to state with any conhdence whether or not 
the smooth lattice method will be of any use in numerical relativity. However the results are 
very encouraging. Furthermore the strong theoretical basis for the method, its central use of 
the equivalence principle and its computational simplicity gives us strong reasons to believe 
that the method will prove to be useful. Of course much more work needs to be done. 

Currently we are adapting these ideas to a 3-1-1 formulation in which the smooth lattice 
method is applied to the 3-geometry while retaining a continuous time coordinate. We will 
employ a dynamical 3-dimensional lattice whose evolution will be controlled by the standard 
ADM 3-1-1 equations (which will be applied directly to the 3-lattice). The smooth lattice 
equations will be used in each time evolution step to aid in the calculation of the source 
terms. For example, at each time step we would be given a discrete 3-metric on a 3-lattice 
from which we could estimate the 3-curvatures by way of the smooth lattice equations. There 
is no role for the Bianchi identities in this phase of the integration, we simply have to live 
with whatever estimates we get for the curvatures. However we can expect the Bianchi 
identities to play an important role in the solution of the initial value problem where the 
initial 3-metric is being constructed. We can speculate that the constraints associated with 
the Bianchi identities on the 3-metric will be preserved by the 3-1-1 evolution equations. This 
is pure speculation at the moment and may in fact prove to be wrong. We will be testing 
this 3-1-1 method on the evolution of the Schwarzschild initial data presented here and we 
hope to report on this work in the near future. 

Another interesting question is how much mismatch in the metric do we incur (or how 
much can we allow) between pairs of adjacent RNC frames? Recall that each RNC frame 
is constructed by choosing its metric such that it agrees with a hnite set of geodesic leg 
lengths. Consider now two adjacent frames that share a common cell. For a typical curve 
in that cell this pair of frames will not yield the same estimates for the length of that curve 
(except on the geodesic edges of that cell). By how much do these estimates differ and is 
there an optimal choice of curvatures that minimises the difference? If such a minimising 
relation exists it must constrain the choice of curvatures in each frame. However, we already 



know that the curvatures must be constrained by the Bianchi identities. How then are these 
constraints related? This too is currently being investigated. 

In both of our examples we chose to solve the Hamiltonian constraint by a radial integration. 
However the constraints are known to be elliptic equations and so should be properly solved 
as a boundary value problem. This leads us to an interesting problem. How can asymptotic 
conditions such as the conformal factor behaves like 1 + k/r for large r be imposed using 
local Riemann normal coordinates? The difficulty is that this asymptotic condition is a 
statement about global properties of the space whereas the Riemann normal coordinates are 
purely local. Another way of looking at this is to ask how can two observers, Joe and Mary, 
each in their own RNC frame with Mary centred on the throat and Joe in the distant flat 
regions, determine which observer is in the asymptoticly flat part of the space. Recall that 
each observer is in a locally freely falling frame and the only geometric quantities that they 
can measure are the curvature components. Thus the statement that observer Joe is in the 
asymptoticly flat region must say something about the curvatures measured by Joe relative 
to those measured by Mary. It would seem that the asymptotic boundary conditions would 
need to be formulated in terms of the curvatures. Yet such a conclusion seems at odds with 
the standard theory of elliptic pde’s. We seem to be imposing boundary conditions on the 
second derivatives of a function which is a solution of a second order elliptic equation. This 
is indeed odd and needs to be resolved if this method is to be applied to general initial value 
problems. 


7. Appendix 

This appendix, which is part of a larger report [5] on Riemann normal coordinates, is 
included here to provide a self-contained derivation of the two smooth lattice equations (1.2) 
and (1.3). This will require some preliminary discussion on Riemann normal coordinates. In 
deriving our equations we will make repeated use of series expansions in the curvature (and 
its derivatives) around flat space. We have already seen this in action when we wrote the 
metric as 

9fip{x) = - ^Rfj,aup -h O(e^) (1.1) 

We shall prove this in section (7.3) below. 

This expansion can be interpreted as a flat space part plus a (presumably) small quadratic 
perturbation. Naturally one might have reservations about the convergence of such an ex¬ 
pansion. We shall first establish that for sufficiently small regions in which this RNC is 
dehned the truncation errors are negligible. We will do this by introducing a conformal 
transformation of the original metric. This will prove to be extremely useful. All of the 
results in this appendix will subsequently be quoted in terms of these conformal coordinates 
whereas in the body of this paper our results are expressed in the physical coordinates. 
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For the remainder of this appendix we will be concerned with one RNC frame. We will 
also refer to the region over which the RNC is dehned as the patch (which in the body of 
the paper consisted of two consecutive cells of the tube or the set of triangles attached to a 
vertex). 

7.1. Conformal coordinates 

Let the typical length scale of the patch containing O be e. Let the coordinates of the patch 
be and let the coordinates of O be x^. Now dehne a new set of coordinates by 

x^ = x^ ^- ey^ 


Then 




Dehne the conformal metric ds by 


gfiv{x)dx^dx'^ 

+ ey)dy^dy’' 


ds^ = g^v{Xi, + ey)dy^dy'' 

= gfiviV: e)dy^dy'' 

From the above it is easy to see that, at O, 

g^iv = g^iv 

gfiu,a ~ ^g^jLV,a / 

g^UjaP ~ ^ g^v,al3 J 

where the partial derivatives on the left are with respect to y and those on the right are with 

respect to x. For each higher derivative an extra power of e will appear. 

From the above we immediately obtain 


R 




e^R 


UvaP 


and as R^pap is independent of e, we see that 


R^luaf3 = 0{e^) 


for e << 1. 

Clearly, then, as e —0 the conformal metric is hat. 

There are now two ways to look at the patch. We can view it (in the original coordinates 
as a patch of length scale e with a curvature independent of e. Or we can view it (in 
the conformal coordinates y^) as a patch of hxed size but with a curvature that varies as e^. 
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This later view is useful since in using it we can be sure that the series expansions around 
flat space are convergent (for a sufficiently small e). 

We will use these conformal coordinates for the remainder of this Appendix. We will also drop 
the tilde and revert to as the generic coordinates (even while working in the conformal 
frame.) 


7.2. Riemann normal coordinates 

In Riemann normal coordinates the geodesics through O must all be of the form 


x^{s) = a^s 


for some set of numbers a^. By direct substitution into the geodesic equation 


0 


d?x^ „ dx^ dx^ 


and its derivatives, one obtains, at the origin O, 


0 

0 


-pfi 

a/3 

pP I pM pM 

^ af3,v ' I3v,a ' ^ va^P 


(7.2.1) 


(7.2.2) 

(7.2.3) 


It is easy to see, by continuing in this way, that all symmetric derivatives of the connection 
vanish at the origin in Riemann normal coordinates. 


7.3. Metric 

Consider a Taylor series expansion of the metric around the origin O, namely, 

+ g^y^ap 2 ^ ) 

There is no linear term because g^y^a = 0 at the origin. It is a simple algebraic exercise to 
show, given (7.2.2) and (7.2.3), that 



^aP,^ ~ 3 {^^aPv + Pav) 

(7.3,1) 

from which it follows that 

and hnally 

9jj,v,aP 2 {R^avP T R^Pva) 

R^vaP ~ 9av,^iP 9a^,vP 

(7.3.2) 
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Substituting these into the above we obtain 


9iJiv{,x) 



+ 0{e^) 


( 1 . 1 ) 


7.4. Connection 

We can also propose a Taylor series expansion for the connection about the origin O, namely, 

Our hrst observations are that 


P rpT 


X^X 
',pr~ 


+ ■■■ 


_ ml 

a(3,p 


= 0{e^ 

= 0(P) 

and in general the n-th derivative of T will be at O. This follows by simple inspection 

of the standard formula for computing the metric connection and the previously stated 
asymptotic behaviour of the conformal metric (7.1.1). 

We are only interested in the leading term in the above expansion, and so after using (7.3.1) 
we obtain 


+ 0(P) 


(7.4.1) 


7.5. Geodesics 

For the geodesics, we employ a series expansion in s, the distance measured along the 
geodesic, 

a;^(s) = Oq + Oj's + 02^ + H- (7.5.1) 

We will defer for the moment stating the nature of the truncation error. Our primary aim 
here is to determine as many of the af as we can in terms of just and R^uap- This can be 
done by demanding that the above expansion for x^{s) is a solution of the geodesic equation 

(7,2,1). 

The basic steps are to substitute (7.5.1) into (7.4.1) and to then substitute all of these 
quantities into the geodesic equation (7.2.1). The result is a polynomial in s and as this 
must be identically zero for all s, we equate the separate coefficients of powers of s to zero. 
For the first two terms and we obtain, respectively, 

0 = «2 + 

0 = Og + ^al3,p + 2000^02^ + ^ai3,pT + 2aQaQa2af j 


- 32 - 



The term is zero in view of (7.2.3). Thus, to order e^, we obtain 

“2 = + 0(€") 

Clearly this process can be developed in full to obtain recurrence relations amongst all of 
the remaining a^. We will not need these but what we do require is their behaviour in e. It 
is not hard to see that in the generic equation for the leading terms will be 

0 = «n + + ■ ■ ■ 

Since we have already established that is 0{e^) it follows that an will be 0{e'^) for n >2. 
Assembling the above results leads dually to 



Notice that Oq and remain undetermined ~ they can only be computed from appropriate 
boundary or initial conditions. 

7.6. Geodesic boundary value problem 

In this case we are looking for the geodesic which passes through two given points. Let the 
coordinates of initial point be and those for the final point be x^. Suppose the geodesic 
distance between the two points is Lij. The cannot be freely specified as they must be 
derivable from the metric and the coordinates. A equation for will be given in a later 
section (7.7). 

Our aim is to solve for Oq and such that 

x'^{s = 0) = xf = Uq + 0{€^) 

x^{s = Lij) = = a[( + a^'Loi + ^R^aPp + 0{e^) 

The first equation is easy to solve, namely, Qq = x^ + 0{e^). However, the second equation 
does appear to pose a bit of a problem - it looks like a nasty quadratic equation for each 
of the a^. Fortunately this equation can be solved by an iterative method to within 0{€^). 
The starting point is to first substitute for Oq to obtain 

+ a^Loi + —R^aPp -^01 + 

We can then re-write this as 

^1=^ - ^R^app a?^01^ + C»(e^) 
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which we will evaluate as a fixed point iteration scheme. 
Since R^app = we obtain the first approximation 




+ 0(e2) 


This can now be substituted back into the previous equation leading to the second approxi¬ 
mation 

ai = ^ - ^R'^app + C>(e^) (7.6.1) 

where Axfj = Xj — xf. This is as far as we can proceed with this iteration scheme because 
the errors in this approximation for ai and those in the equation we are iterating on are 
both 0{e^). Combining these results for Oq and and substituting into (7.5.2) leads to the 
following equation for the geodesic passing through the two points x^ and x^ 


x^is) = xf + \/\x% - <AxgAxg. + 0(i) (7.6.2) 


where A = s/Lij. 


7.7. Geodesic distance 


Consider two points with coordinates x^ and x^. Since there exists, by assumption, a unique 
geodesic joining this pair of points, the distance between them should also be uniquely defined 
in terms of their coordinates and the metric. 

Our aim is to evaluate, along the geodesic. 




dx^ dx^ A 

d\ d\ J 


1/2 


dX 


The equation for x^{\) is simply (7.6.2) for 0 < A < 1. This can be substituted into the 
expansion (1.1) for g^v{x) with the result 


gp,y{x{X)) = gf,u - ^Rpai^p {xf + XAxfj) + AAxg.) + C>(e^) 


It is a simple matter to substitute these into the integrand, leading to 




l^RpauP xfx^AxljAx'^ij + 0{ 


The important point to note is that this result does not depend on A. Thus the integrand is 
constant and so the integration is trivial. The result follows immediately, 

= g^i^AxfjAx^j - ^R^ai^g xfx^Ax^^-Ax^j + 0{e^) (7.7.1) 
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From this result it is easy to establish, using the symmetries of R^pap-, fhe following equivalent 
equations for L‘f- 

L% = xfj)Ax^jAx'^j + 0{e^) 

= 9tiuAxfjAxij - ^R^aup xfx'^x'^x'^j + 0{e^) 

= 9^^uAx^JAx’^j - ^R^aup Ax^iAx1^Ax>^-Axlj + O(e^) (7.7.2) 

where x^j = [xf + Xj )/2 and where Xq are the coordinates of the origin (which in onr case 
are zero). 

Note that when these eqnations are re-expressed in terms of the physical metric all that 
changes is that the error term C7(e^) is replaced with O(e^). 

7.8. Generalised Cosine law 

Consider a geodesic triangle with vertices i,j and k. We would like to be able to compute the 
angles subtended at each vertex in terms of the usnal quantities, the metric, the coordinates 
etc. We will develop the appropriate eqnations in two stages. First, we will consider the 
simple case of compnting the angle at a vertex coincident with the origin. We shall then 
generalise this result to the case were all three vertices are distinct from the origin. 

To start the ball rolling consider a geodesic triangle with vertices i,j and O, the origin. We 
seek an eqnation for the angle between the geodesic segments joining O to i and O to j. The 
nnit tangent vectors to these geodesic segments at O are, from eqnation (7.6.2), 

^oi = ^^oi/Loi 
^oj = ^^oj/Loj 

Now let 9 ij be the angle snbtended at O. Then 

cos^y = = 9^^uAx^^Ax’'oj/{LoiLoj) 

We can obtain two usefnl variants of this equation by writing, first, AxC = AxC — AxC and 
second, AxC = AxC -|- Axfj. This gives 


LoiLoj cos 9ij = g^p ~ 

= dfiuAx^^ (Ax^j + Ax^j) 

Adding these two equations leads to 

2LoiLoj cos 9ij = Lli + Llj - g^pAx^-Ax'lj (7.8.1) 
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However, from equation (7.7.2) we see that 


= Ljj + Ax'^^Ax'^^iAx'^jAx^j + 0{e^) 

Thus we have 

2LoiLoj cos Oij = Lli + Llj - Ax" + 0{e^) (7.8.2) 

With this equation we have achieved our hrst aim : to obtain an equation when the vertex 
resides at the origin. To obtain an equation applicable to the general case, where the vertex is 
not at the origin, we can imagine transforming to a second set of Riemann normal coordinates 
with an origin at some other point, say O'. We can do this simply by shifting the coordinates, 
eg. x^ ^ x^+c^. The coordinates, metric and Riemann components at the respective origins 
will therefore be related by 


x'^ = x^ + 

g't,AO’) = 9^.(0) - cV + C>(e3) 

+ 0(X) 

Now the important observation is that the above eqnation (7.8.2) is covariant with respect 
to this transformation (whereas (7.8.1) is not). That is, it applies to any three vertices of a 
geodesic triangle. Let us now relabel the vertices as i,j and k. Then the angle snbtended at 
vertex k can be computed from 

‘2LikLjf; cosOij = — -R^a^g Ax^^Ax^j^Axjj^Ax^-j, + 0{e ) (7.8.3) 

As with the geodesic length eqnation (7.7.1), one need only replace the 0{e^) error term 
with 0{€^) to obtain the correct equation in the original physical coordinates. 
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Figure 1 A simple approximation to a 2- 
sphere. Though the legs of the tetrahedron 
appear straight they are taken to be geodesic 
segments of the 2-sphere. 


Figure 2 Better approximations for the 2- 
sphere are obtained by sub-dividing each tri¬ 
angle according to this pattern. Newly cre¬ 
ated vertices are placed at the centre of old 
legs and are later displaced radially so that 
they, the new vertices, touch the 2-sphere. 



re error m rj as a runction oi the 


R. Note that the errors vanish a 

■g I- 










1 + ri/{rp^) 



Figure 6 This demonstrates that even though th 
unbounded for large /, they do vanish as i? —2 at 





edges are global geodesics. 





Figure 8 The generic pair of cells for one Riemann normal coordinate frame. The ; 2 -axis runs 
up the middle of the tube. 





Figure 10 Cross sections of Figure 9 for various values of 1. In each case the error 
converges quadraticly with respect to e. The irregular behaviour in ER^ is probably 
due to round-off error. The curves correspond to / = 5,10, 20,40 and 80 with the 
error generally increasing with 1. This applies also to Figure 12. 





or th( 


"hich the a 
n ERx are 
y small. 






;ure 12 Cross sections of Figure 11. In each case the error converges quad 
ih respect to e. The irregular behaviour in ER^ is due to the oscillations 
•vatures. 














